'''
Created on Jun 15, 2010

@author: oabalbin
'''

import sys
import copy
import pickle
from optparse import OptionParser
from datetime import datetime

import amarillo.copa.copa_analysis as ca
import amarillo.copa.subset_genes as sg

if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", "--outfolder", dest="outfolder",
                            help="output folder")
    optionparser.add_option("-s", "--normalsFile", dest="normalsFile",
                            help="annotation file with the normal samples")
    optionparser.add_option("-p", action="store_true", dest="dumpfile",
                            help="use a dump file of the expression matrix. It is passed in annotFile")
    optionparser.add_option("-g", "--shortgenelist", dest="shortgenelist",
                            help="annotation file with a short genelist")
    optionparser.add_option("-t", "--tissueName", dest="tissueName",
                            help="name of the tissue analized")


 

    
    (options, args) = optionparser.parse_args()
    
    
    t = datetime.now()
    tstamp = t.strftime("%Y_%m_%d_%H_%M")
    
    outfiles = ['outfile_outliers','outfile_outliers_names','outfile_outliers_list','outfile_outliers_matrix','outfile_copa_matrix',
                'outfile_foldchange_matrix','outfile_foldchangesinsample']
    time_stamp_outfiles=[]
    for name in outfiles: 
        time_stamp_outfiles.append( options.outfolder+tstamp+'_'+options.tissueName+'_'+name+'.txt')
        
    thisfile = open(options.annotFile)
    outfile_outliers= open(time_stamp_outfiles[0],'w')
    outfile_outliersinsample= open(time_stamp_outfiles[1],'w')
    outfile_outliers_list = open(time_stamp_outfiles[2],'w')
    outfile_outliers_matrix = open(time_stamp_outfiles[3],'w')
    array_dump_file = options.annotFile + '_dump_file'
    outfile_copa_matrix = open(time_stamp_outfiles[4],'w')
    outfile_foldchange_matrix = open(time_stamp_outfiles[5],'w')
    outfile_foldchangesinsample = open(time_stamp_outfiles[6],'w')
    ####################
    ##### Parameters
    
    normal_samples_list = ca.list_of_names(open(options.normalsFile))
    
    #percentils_list = [50]
    percentils_list = [75,80,85,90] #[80,85,90,95,98] #[75,80,85,90] #[80,85,90,95,98] 
    method='mean' # median, mean, sum
    outliers_cutoff=200 #100 #200
    nout2report=20
    expression_cutoff=0.0       # minimum median RPKM of expression
    base_expression_cutoff= 50  # In percentil. It is used to report the expression value of a particular gen in a sample. 
                                # The median usually can be used, however because in many cases median=1, It is maybe desireable to use for e.g the75%. 
    median_shift=1.0
    sdv_factor=0.25             #0.25
    not_report_isoforms=False   # Report gene outlier isoforms
    
    ########### Program
    
    if options.dumpfile:
        thisGeneArray = pickle.load(open(options.annotFile))
    else:    
        thisGeneArray = ca.create_myArray(thisfile, expression_cutoff, normal_samples_list)
        pickle.dump(thisGeneArray,open(array_dump_file,'w'))
        
    if options.shortgenelist:
        allGeneArray = copy.copy(thisGeneArray)
        new_genelist = ca.dict_of_names(open(options.shortgenelist))
        thisGeneArray = sg.create_array4subset_genes(thisGeneArray, new_genelist)
        print thisGeneArray.expVal.shape

    
    
    thisGeneArray.copa_norm_across_samples(median_shift)
    thisGeneArray.copa_norm_in_sample()
    thisGeneArray.fold_change_accross_samples()
    #outliers_cutoff = thisGeneArray.expVal.shape[0]
    #print outliers_cutoff

    ############# Print expression matrix after Copa normalization
    ca.print_copa_matrix(thisGeneArray, outfile_copa_matrix)
    #ca.print_foldchange_matrix(thisGeneArray, outfile_foldchange_matrix)
    #ca.print_foldchange_outliers_by_samples(thisGeneArray,outfile_foldchangesinsample, nout2report)
    #sys.exit(0)
    ############ Run the outlier analysis for different percentiles
    
    gene_outliersample_dict, copa_rank_dict, copa_score_dict = ca.copa_at_different_percentiles(thisGeneArray,percentils_list)

    ############ Summarize the results from the different percentiles
    
    final_list_genes, gene_copascore_values = ca.summarize_results_copa_percentiles(thisGeneArray,
                                                                                  gene_outliersample_dict, copa_rank_dict, 
                                                                                  copa_score_dict, percentils_list, sdv_factor, method)
    
    ############# Identify all outliers for each sample
    
    list_outliers_report, outliers_in_samples, outliers_in_samples_copascore = ca.identify_outliers_in_sample(thisGeneArray,final_list_genes, gene_outliersample_dict, 
                                                                                                              gene_copascore_values, outliers_cutoff, not_report_isoforms)
    
    
    ### Print outlier list and outlier matrix. 
    ca.print_list_of_outliers(list_outliers_report, outfile_outliers_list)
    ca.print_matrix_of_outliers(thisGeneArray.get_list_of_samples(), list_outliers_report, outliers_in_samples_copascore, outfile_outliers_matrix)
            
    
    ############# Print sample : outliers gene in that sample
    
    ca.print_outliers_by_samples(thisGeneArray, gene_copascore_values, outliers_in_samples, outfile_outliersinsample, nout2report)
    